Controlling the position of traveling waves in reaction-diffusion systems 
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We present a method to control the position as a function of time of one-dimensional traveling 
wave solutions to reaction-diffusion systems according to a pre-specified protocol of movement. 
Given this protocol, the control function is found as the solution of a perturbatively derived integral 
equation. Two cases are considered. First, we derive an analytical expression for the space (x) and 
time (t) dependent control function / (x,t) that is valid for arbitrary protocols and many reaction- 
diffusion systems. Second, for stationary control of traveling waves in one-component systems, the 
integral equation reduces to a Fredholm integral equation of the first kind. In both cases, the control 
can be expressed in terms of the uncontrolled wave profile and its propagation velocity, rendering 
detailed knowledge of the reaction kinetics unnecessary. In two spatial dimensions, an extension of 
the proposed method allows to control the shape of traveling waves. 



Introduction. Many control methods were devised for 
reaction-diffusion systems (RDS), e.g. global periodic 
forcing, feedback controls, time delayed feedback etc. 
[l|| . For example, unstable patterns can be stabilized 
by global feedback control, as was shown in experiments 
with the light-sensitive Belousov-Zhabotinsky reaction 
[2|. Varying the global light intensity forces a spiral 
wave tip to describe a wide range of open and closed 
hypocycloidal trajectories [3J. A spatiotemporal feedback 
control guiding a propagating wave along a desired tra- 
jectory was realized experimentally in |J|. Dragging of 
a traveling chemical pulse [5| on an adressable catalyst 
surface [6j was accomplished by a moving, localized tem- 
perature heterogeneity. 

We consider perturbed reaction-diffusion systems of the 
form 



d t u = Dd^u + R (u) + eG (u) f (x, t) , 



(1) 



where D is a diagonal matrix of constant diffusion co- 
efficients, f is a spatiotemporal perturbation coupling 
via a possibly u-dependent symmetric matrix Q and R 
is a typically nonlinear reaction function. The unper- 
turbed (e = 0) solution is assumed to be a traveling wave 
U c (?) , ? = x — ct, stationary in the reference frame co- 
moving with velocity c, so that 



DdpJc (?) 



c%U c (?)+R(U c (?)) = 0. 



The eigenvalues of the linear operator 



C = Ddl + cd/: 



2?R(U C (?)) 



(2) 



(3) 



determine the stability of the traveling wave, with 
PR (U c (?) ) the Jacobi matrix of R evaluated at U c (?) . 
We always assume U c to be stable so that A = is the 
eigenvalue with largest real part and the Goldstone mode 
W (?) = TJ' C (?) the corresponding eigenfunction. Fur- 
thermore, we assume that C has a spectral gap so that 
the eigenvalue with next largest real part is a finite dis- 
tance away from the zero eigenvalue. Because C is in 



general not selfadjoint with respect to the standard in- 
ner product in function space, the eigenfunction W^ (£) 
of the adjoint operator & to eigenvalue zero, the so-called 
response function, is not identical to W (?). 
Perturbatively for < e < 1, an equation of motion for 
the position cp (t) of the traveling wave under perturba- 
tions can be derived, 



4> (*) = c - 



K c 



W* T (x) G (U c (x)) f(x + <f>(t), t) da:, 

(4) 



with initial condition (j) (t ) — 4>$ and 



K C = (W\\J' C )= [ W tT (x) U' c (x) dx. 



(5) 



For a derivation see e.g. [7( and [7|, [8( for applications 
of Eq. ([4]). Similar equations of motion were derived 
for spiral waves in two spatial dimensions [9||. A related 
approach is the direct soliton perturbation theory [lOj 
developed for conservative systems supporting traveling 
waves as e.g. the Korteweg-de Vries equation. 
Without exception, we assume e = 1, meaning in 
turn that some norm of f, as e.g. the uniform norm 
| |f (cc, t)\ 1^ = sup^gjj |f (x,i)|, must be small for all times 
for Eq. (j4|) to be valid. For monotonously decreasing 
front solutions, we define its position as the point of 
largest slope, while for pulse solutions it is the point of 
maximum amplitude of an arbitrary component. 
In this paper, we do not consider Eq. (j4]) as an ODE for 
the position of the wave under the perturbation f , but 
as an integral equation for the control f realizing an arbi- 
trary but given protocol of movement <fi (t) . We assume 
that the wave moves unperturbed until reaching position 
</>o at time to, upon which the control is switched on. 

Control realizing an arbitrary protocol of movement. 
A general solution for the control function f with arbi- 




We assume that the concentrations Ci/ 2 of the chemi- 
cal species A 1 / 2 , Ci/ 2 = [A 1 / 2 ], can be controlled spa- 
tiotemporally, C\ i 2 — ¥ C\i 2 + ef (x, t), giving an additive 
control function (Q (u) = 1) in the case of control via 
C 2 and a multiplicative control (Q (u) — u 2 ) in the case 
of control via C\. In the bistable parameter regime, the 
unperturbed one-dimensional Schlogl model has an ana- 
lytically known traveling front solution U c [13fl connecting 
the stable upper and lower homogeneous steady states as 
x — > ±00. 

We want to move the front back and forth in a sinusoidal 
manner with amplitude A and period T, 



Figure 1. Multiplicative position control of a Schlogl front 
solution. The control function couples to the quadratic term 
of the Schlogl model. Comparison of numerically obtained 
position over time data (red dashed line) and protocol (black 
solid line). Inset shows velocity. See also SI in |ll| . 



trary function h (x) and arbitrary protocol <f> (t) is 

f (as, t)=(e-4> (*)) ^Q- 1 (U c {x-4> («))) h (x - 4> (t)) ■ 



G, 



W tT (x) h (x) dx, 



(G) 



(7) 



where Q" 1 denotes the matrix inverse to Q which is as- 
sumed to exist. A control proportional to the Goldstone 
mode prevents large deformations of the wave profile 
and shifts the wave profile as a whole ll|. We choose 
h (x) = V' c (x) and 



f (x, t)=(c-j> (t)) G~ x (U c (x-<f> (*))) U^ (x - 



4>{t)). 

(8) 



The analytical protocol is compared with position over 
time data obtained by numerical simulations of the con- 
trolled RDS with no flux or periodic boundary conditions. 
The uncontrolled traveling wave is used as the initial con- 
dition. A large amplitude of U^. renders the equation of 
motion Eq. ([4]) invalid, but can be counteracted by a 
choice of protocols with a velocity <\> close to the velocity 
c of the uncontrolled traveling wave. 

Multiplicative control of a Schlogl front. The Schlogl 
model was proposed in [T2| as a chemical reaction mech- 
anism 



A x 



2X ^ 3A, 

K 



X%A 2 



(9) 



It gives a cubic reaction function for the time evolution 



of the concentration u of X, u — [X], (with kf, 2 = 1) 



R(u) =CiU 2 -u 3 -u + C 2 - 



(10) 



(j) (t) = A + A sin 



2tt 



t + Ai 



(11) 



Aq and 



via a spatiotemporal control of C\ {Q (u) = 

A\ are determined by <f> (to) = 4>o, <fi' (to) = c so that the 

protocol is continously differentiable at t = to- 

The position of the numerical front solution follows the 

protocol very closely, see Fig. [T] Note that the maximum 

enforced front velocity, maxt (/)' (t) = 7.854, is many times 

larger than the velocity c = 0.660 of the uncontrolled 

front, see inset of Fig. [TJ 

FitzHugh-Nagumo model. We apply the position con- 
trol to the stable traveling pulse solution of the FitzHugh- 
Nagumo model [l4| 



dtu =D u d x u + 3u — u — v + ef u (x, t) . 
d t v =D v d 2 x v + e(u-S)+ ef v (x, t) . 



(12) 



We assume to have the ability to control two independent 
additive parameters spatiotemporally, one parameter in 
each the equations for the activator u and inhibitor v so 
that Q (u) — 1. As an example, we consider a protocol 
4> (t) which reverses the propagation direction of the un- 
perturbed pulse in Fig. [2] (see also movie S2 in [ll|). The 
protocol is constructed smoothly of linear and sinusoidal 
parts and periodic boundary conditions are applied. Pro- 
tocols with a velocity jump at t = to can lead to unde- 
sired behavior. For example, a protocol which tries to 
stop the pulse instantaneously can lead to a stationary 
pulse emitting traveling waves to the right, as is shown 



in the supplemental material S3 in [11 1 



Shaping a traveling wave. For a two-dimensional ex- 
tension of the equation of motion, we allow for devia- 
tions in the lateral direction, <f> — <f>(t,y). By a multi- 
ple scale perturbation expansion around the unperturbed 
plane wave U c traveling in the x-direction, we derive 
the perturbed nonlinear phase diffusion equation |15j 
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Figure 2. Reversal of propagation direction of a FitzHugh- 
Nagumo pulse. The position of the controlled pulse (red 
dashed line) follows the protocol (black solid line) very closely. 
Insets show snapshots of the controlled pulse profiles (left) of 
activator (blue solid line) and inhibitor (purple dotted line) 
and the corresponding control functions f u and f v (right). 
See also movie S2 in |T1| . 



c+^ + ad," 



K, 



W^ T {x)g{\J c {x))f{x + <p,t)Ax (13) 



1 



with a = --— (Wt,DUj.). No flux boundary conditions 

in the y-direction of the two-dimensional RDS correspond 
to Neumann boundary conditions for 4>. A solution for 
the control function is 
c 

2< 

(14) 



f(s)=(c 



+ -do 12 + ad," g- 1 (U c (x - <j>)) U' c (x 



which we apply via an additive control scheme (Q (u) — 
1) to the Schlogl front solution where a = D. See Fig. 
[3] for a snapshot and S5 in [11] for a movie of the time 
evolution of a front under a control enforcing the tran- 
sition from an unperturbed plane wave traveling in the 
x-direction to a sinusoidally shaped stationary front. 

Stationary control of traveling waves. For one- 
dimensional traveling wave solutions of one-component 
systems, the general expression for the adjoint Goldstone 
mode is 



W^0=e ci/D K(0- 



(15) 



Introducing the inverse function T = <f> , one can for- 
mulate a Fredholm integral equation of the first kind for 



fl(0)= c 



T'(0) 



Kr 



K\ 



x)f(x)dx, (16) 



Figure 3. Shaping a Schlogl front according to a protocol 
4>(t,y) (black line). Yellow (red) corresponds to a high (low) 
value of u. See also movie S5 in Hill. 



with kernel K (x) = &- cx l D U' c (-x) Q {U c (-x)) and in- 
homogeneity g. Eq. (JTHJ) can be solved with the help of 
the convolution theorem for the two-sided Laplace trans- 
form. 

As an example, we choose a protocol which drives the 
propagation velocity c > to zero, 

4> (*) = | (1 + tanh (k (ti - t))) ,ti>t 0f k> 0. (17) 

In the limit k — > oo, we recover a protocol which instan- 
taneously stops the wave at t = t\, 



lim d? (t) = cO (ti - t) , 

k— »-oo 



(18) 



where represents the Heaviside Theta function. The 
inhomogeneity g is determined as 



' (0) = K c 



exp(^(ci o + 0-</> o )) 



,2H„ 



.,2fc«i 



(19) 



The coupling function is set to Q (u) = 1. 
As the model, we choose the Schlogl model Eq. ([TO)) in 
the rescaled form R (u) ~ —u (u — a) (u — 1) with front 
solution U c (£) = 1/ (l + CX P v»/V*)) anc l velocity c = 

—= (1 — 2a) [13j for D = 1. The region of convergence of 

v 2 

the Laplace transforms of kernel K and inhomogeneity g 

determine the allowed range of values for k, 



< k< C -(^=-c 
2 \V2 



(20) 



This amounts to a minimum acceleration (or maximum 
deceleration) at time t = t\, 



o<0(M = ^<-(|) 2 (^- C 



(21) 
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Figure 4. Deceleration of a Schlogl front by an additive sta- 
tionary control. Red dashed line is the result of numerical 
simulations, black solid line is the protocol. Shown is the po- 
sition <j> and velocity 4>' (right inset). The front profile (blue 
solid line) is slightly deformed in the region where the control 
(purple dotted line) is large, see left inset. See also S9 in [ill ). 



which can be realized for such a protocol. For the control 
function, we find 



K r c sin 



/(z) = - 



VW C 2 



2k 



exp 



■2k 



O-0o)) 



^(l+e^i-io)) ( c 2 + 2fc) 



(22) 



The control is diverging to — oo as x — > oo, which we cir- 
cumvent by clipping / so that locally R (u) + ef (x) = 
has three real roots and bistability is preserved at every 
point in space. A more systematic approach to prevent 
a diverging solution would be to consider the Fredholm 
integral equation Eq. (| 16[) supplemented with inequality 
constraints for the control function. 
Under the control Eq. (|22|) , the velocity of the numerical 
solution first follows the velocity protocol closely. Devia- 
tions arise when the transition region of the front enters 
the domain where the absolute value of the control is 
large, see right inset of Fig. [5] This difference in the 
velocities leads to an accumulated difference for the po- 
sition at which the front is stopped. The numerical front 
profile is slightly deformed in the region where the control 
is large, see left inset in Fig. @] 

Conclusions. The proposed control methods work for 
a large class of systems supporting traveling wave so- 
lutions and many different protocols. To construct the 
control functions, the profile U c of the uncontrolled trav- 
eling wave must be known sufficiently accurate to deter- 
mine its derivative U^,. Furthermore, knowledge of the 
velocity c as well as the coupling matrix Q (u) is nec- 
essary. The value of the diffusion coefficient D is re- 
quired for stationary control, Eq. (|16p , and spatiotempo- 
ral two-dimensional control, Eq. (|14|l. of traveling waves 
in one-component systems. In these cases, knowledge 
of the nonlinearity R (u) is not necessary, which ren- 



ders the control methods useful for applications where 
models are only approximately known. Spatiotemporal 
two-dimensional control of multi-component systems re- 
quires the response function W^ to determine a. It is not 
known how to express W^ in terms of the uncontrolled 
wave profile XJ' C and the velocity c for the general multi- 
component case. 

The one- and two-dimensional spatiotemporal control so- 
lutions, Eq. ([8j) and Eq. (fT4|) . shift the wave profile ac- 
cording to the protocol and deform the wave profiles only 
weakly. Because traveling wave profiles U c (£) usually 
decay exponentially fast as £ — > ±oo, the spatiotemporal 
controls in one and two spatial dimensions, Eq. ([8]) and 
Eq. (fT4| . are localized in the x-direction. 
The solution for the one-dimensional spatiotemporal con- 
trol, Eq. ©, can also be seen as a so-called singular 
solution to an u nreg ularized and unconstrained optimal 
control problem [16j : solve the RDS Eq. (JlJ for the con- 
trol f , substitute u (x, t) — U c (x — <fi (t)) as the desired 
distribution and use the property of traveling wave so- 
lutions Eq. ([2|) to eliminate the diffusion and reaction 
terms. 

If the control affects only the equation for, say, the in- 
hibitor in the FitzHugh-Nagumo model Eq. (|T2j) . then 




S{u) = 



1 



is not invertible. Note that one can 



still find a general solution to the integral equation with 
arbitrary function h as (with WJ being the inhibitor com- 
ponent of the adjoint Goldstone mode) 



fv (X, t) 



G, : 



"i/.J^M*- 



Wl (x) h (x) dx. 



0(t)), (23) 



(24) 



However, the choice of h (x) = V' c (x) for the FitzHugh- 
Nagumo model leads to the undesired spontaneous gener- 
ation of pulses, see S4 in [llj. The assumed invertibility 
of Q excludes control schemes with a number of inde- 
pendent control parameters smaller than the number of 
components of u. 

Reliability of the proposed control methods is difficult 
to assure. A necessary condition on the traveling wave 
solution is the existence of a spectral gap. For example, 
applying a spatiotemporal control to the marginally sta- 
ble front solution of the Fisher equation [17| , which does 
not exhibit a spectral gap [18], leads to a front profile 



growing indefinitely to — oo, see S6 and S7 in [llj . Large 
control functions can destroy the traveling wave and of- 
ten lead to the spontaneous generation of pulses due to 
the underlying excitability of the system, as was also ob- 
served in [5|. The three component Oregonator model 
[19j . for example, appears to be very sensitive, allowing 
only protocols with velocity </>' (£) close to the velocity c 
of the uncontrolled wave, see S8 in [ll|. It is desirable 
to establish estimates how large a control is allowed to 



be so that the wave's position is controlled in a reliable 
way. Protocols should be smooth, although protocols 
with a discontinuous velocity work often in the case of 
the Schlogl model. 

Further generalizations to two spatial dimensions, other 
systems supporting traveling waves as e.g. the Korteweg- 
de Vries equation and the relation of the proposed control 
methods to optimal control will be studied in a forthcom- 
ing publication. 
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